Figures

Luhrmann, Weisman, et al.

This notebook generates plots for all Studies.

source("../scripts_general/dependencies.R")
source("../scripts_general/custom_funs.R")
source("../scripts_general/var_recode_contrast.R")
source("../scripts_general/data_load.R")

# rescale to 0-1
d1_fig <- d1 %>%
  filter(!is.na(country)) %>%
  mutate(spev_score = spev_score/1,
         pv_score = pv_score/3,
         abs_score = abs_score/1)

d2_fig <- d2 %>%
  filter(!is.na(country)) %>%
  mutate(spev_score = spev_score/1,
         por_score = por_score/1)

d3_fig <- d3 %>%
  filter(!is.na(country)) %>%
  mutate(spev_score = spev_score/4,
         dse_score = dse_score/5,
         abs_score = abs_score/1)

d4_fig <- d4 %>%
  filter(!is.na(country)) %>%
  mutate(spev_score = spev_score/4,
         dse_score = dse_score/5,
         pv_score = pv_score/3,
         por_score = por_score/2,
         abs_score = abs_score/1,
         hall_score = hall_score/3,
         para_score = para_score/1,
         cog_score = (cog_score + 2)/4,
         ctl_score = (ctl_score + 3)/6)

Figure 1 (main text)

Setup

d_all_fig <- d1_fig %>% 
              mutate(subject_id = paste("s1", subject_id, sep = "_")) %>%
  select(study, country, site, religion, subject_id, 
         spev_score, pv_score, abs_score) %>%
  full_join(d2_fig %>%
              mutate(subject_id = paste("s2", subject_id, sep = "_"),
                     site = "urban") %>%
              select(study, country, site, religion, subject_id,
                     spev_score, por_score)) %>%
  full_join(d3_fig %>%
              mutate(subject_id = paste("s3", subject_id, sep = "_"),
                     site = "urban", 
                     religion = "undergraduates") %>%
              select(study, country, site, religion, subject_id,
                     spev_score, abs_score)) %>%
  full_join(d4_fig %>%
              mutate(subject_id = paste("s4", subject_id, sep = "_"),
                     site = "urban", 
                     religion = "undergraduates") %>%
              select(study, country, site, religion, subject_id,
                     spev_score, pv_score, por_score, abs_score)) %>%
  mutate(religion = factor(religion,
                           levels = c("local", "charismatic",
                                      "general population", "undergraduates"),
                           labels = c("Faith of local salience",
                                      "Charis. evang. Christianity",
                                      "General population", "Undergraduates")),
         site = factor(site, 
                       levels = c("urban", "rural"),
                       labels = c("Urban", "Rural")),
         study = gsub("study", "Study", study))

Joining, by = c("study", "country", "site", "religion", "subject_id", "spev_score")
Column `site` joining factor and character vector, coercing into character vectorColumn `religion` joining factors with different levels, coercing to character vectorJoining, by = c("study", "country", "site", "religion", "subject_id", "spev_score", "abs_score")
Joining, by = c("study", "country", "site", "religion", "subject_id", "spev_score", "pv_score", "abs_score", "por_score")

d_all_sum <- d_all_fig %>%
  group_by(study, country, religion, site) %>%
  mutate(total_n = n()) %>%
  group_by(total_n, add = T) %>%
  summarise_at(vars(ends_with("_score")),
               funs(mean = mean(., na.rm = T), 
                    sd = sd(., na.rm = T))) %>%
  ungroup() %>%
  full_join(d_all_fig %>%
              filter(!is.na(spev_score)) %>%
              group_by(study, country, religion, site) %>%
              mutate(spev_score_n = n()) %>%
              ungroup()) %>%
  full_join(d_all_fig %>%
              filter(!is.na(pv_score)) %>%
              group_by(study, country, religion, site) %>%
              mutate(pv_score_n = n()) %>%
              ungroup()) %>%
  full_join(d_all_fig %>%
              filter(!is.na(por_score)) %>%
              group_by(study, country, religion, site) %>%
              mutate(por_score_n = n()) %>%
              ungroup()) %>%
  full_join(d_all_fig %>%
              filter(!is.na(abs_score)) %>%
              group_by(study, country, religion, site) %>%
              mutate(abs_score_n = n()) %>%
              ungroup())

Joining, by = c("study", "country", "religion", "site")
Joining, by = c("study", "country", "religion", "site", "subject_id", "spev_score", "pv_score", "abs_score", "por_score")
Joining, by = c("study", "country", "religion", "site", "subject_id", "spev_score", "pv_score", "abs_score", "por_score")
Joining, by = c("study", "country", "religion", "site", "subject_id", "spev_score", "pv_score", "abs_score", "por_score")

fig1_fun_v2_panel <- function(df = d_all_fig, df_sum = d_all_sum, 
                              which_religion, y, include_text = F,
                              d_width = 1, j_height = 0.02) {
  plot <- df %>% 
    filter(religion == which_religion) %>%
    mutate(facet_lab = interaction(study, site, sep = ": "),
           country = factor(country, levels = levels_country,
                            labels = c("US", "Gh.", "Th.", "Ch.", "Va."))) %>%
    ggplot(aes(x = country, y = !!sym(y), 
               color = country, fill = country)) +
    facet_grid(~ reorder(facet_lab, as.numeric(factor(study)))) +
    geom_point(position = position_jitterdodge(jitter.width = d_width/1.01,
                                               jitter.height = j_height,
                                               dodge.width = d_width), 
               alpha = 0.15, show.legend = F) +
    geom_pointrange(data = . %>%
                      group_by(study, country, religion, site, facet_lab) %>%
                      summarise(mean = mean(!!sym(y), na.rm = T),
                                sd = sd(!!sym(y), na.rm = T)) %>%
                      ungroup(),
                    aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                    position = position_dodge(width = d_width),
                    # shape = 21,
                    fatten = 1,
                    color = "black") +
    scale_color_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2") +
    scale_y_continuous(limits = c(0 - j_height, 1 + j_height), 
                       breaks = seq(0, 1, 0.25)) +
    theme(plot.title = element_text(hjust = 0.5),
          legend.position = "bottom") +
    # theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
    guides(color = F, fill = F, 
           shape = guide_legend(override.aes = list(fill = "black"))) +
    labs(title = gsub("Charis.", "Charismatic",
                      gsub("evang.", "evangelical", which_religion)),
         x = "Country", 
         shape = "Site")
  
  if (include_text) {
    plot <- plot +
      geom_text(data = df_sum %>%
                  filter(religion == which_religion) %>%
                  mutate(facet_lab = interaction(study, site, sep = ": "),
                         country = factor(country, levels = levels_country,
                                          labels = c("US", "Gh.", "Th.", 
                                                     "Ch.", "Va."))) %>%
                  select(study, country, religion, site, facet_lab,
                         !!sym(paste0(y, "_mean")), !!sym(paste0(y, "_sd"))) %>%
                  mutate_at(vars(-study, -country, -religion, -site, -facet_lab), 
                            funs(format(round(., 2), nsmall = 2))) %>%
                  mutate(lab = case_when(
                    grepl("na", tolower(!!sym(paste0(y, "_mean")))) |
                      grepl("na", tolower(!!sym(paste0(y, "_sd")))) ~ "",
                    TRUE ~ paste0(!!sym(paste0(y, "_mean")), 
                                  "\n(", 
                                  !!sym(paste0(y, "_sd")), 
                                  ")"))),
                aes(y = 1, label = lab),
                position = position_dodge(width = d_width),
                color = "black", size = 2.5, vjust = 1)
  }
  
  return(plot)
}

fig1_fun_v2_row <- function(which_df = d_all_fig, which_df_sum = d_all_sum,
                            which_y, ylab, 
                            include_xlab = F, include_title = F) {
  
  p1 <- fig1_fun_v2_panel(df = which_df,
                          df_sum = which_df_sum,
                          which_religion = "Faith of local salience",
                          y = which_y) + labs(y = ylab)
  
  p2 <- fig1_fun_v2_panel(df = which_df,
                          df_sum = which_df_sum,
                          which_religion = "Charis. evang. Christianity",
                          y = which_y) + labs(y = ylab)
  
  p3 <- fig1_fun_v2_panel(df = which_df,
                          df_sum = which_df_sum,
                          which_religion = "General population",
                          y = which_y) + labs(y = ylab)
  
  p4 <- fig1_fun_v2_panel(df = which_df,
                          df_sum = which_df_sum,
                          which_religion = "Undergraduates",
                          y = which_y) + labs(y = ylab)
  
  if (!include_title) {
    p1 <- p1 + labs(title = NULL)
    p2 <- p2 + labs(title = NULL)
    p3 <- p3 + labs(title = NULL)
    p4 <- p4 + labs(title = NULL)
  }
  
  if (!include_xlab) {
    p1 <- p1 + labs(x = NULL)
    p2 <- p2 + labs(x = NULL)
    p3 <- p3 + labs(x = NULL)
    p4 <- p4 + labs(x = NULL)
  }
  
  fig <- plot_grid(p1, 
                   p2 + theme(axis.ticks.y = element_blank(),
                              axis.text.y = element_blank(),
                              axis.title.y = element_blank()), 
                   p3 + theme(axis.ticks.y = element_blank(),
                              axis.text.y = element_blank(),
                              axis.title.y = element_blank()), 
                   p4 + theme(axis.ticks.y = element_blank(),
                              axis.text.y = element_blank(),
                              axis.title.y = element_blank()), 
                   nrow = 1, rel_widths = c(2.3, 3, 1, 2))
  
  return(fig)
}

fig1_row1 <- fig1_fun_v2_row(which_y = "spev_score", ylab = "Spiritual Events", 
                             include_title = F)

Removed 3 rows containing missing values (geom_point).Removed 1 rows containing missing values (geom_point).Removed 1 rows containing missing values (geom_segment).Removed 2 rows containing missing values (geom_point).

fig1_row2 <- fig1_fun_v2_row(which_y = "pv_score", ylab = "Porosity Vignettes")

Removed 11 rows containing missing values (geom_point).no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -InfRemoved 265 rows containing missing values (geom_point).Removed 5 rows containing missing values (geom_pointrange).no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -InfRemoved 766 rows containing missing values (geom_point).Removed 5 rows containing missing values (geom_pointrange).no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -InfRemoved 519 rows containing missing values (geom_point).Removed 5 rows containing missing values (geom_pointrange).

fig1_row3 <- fig1_fun_v2_row(which_y = "por_score", ylab = "Porosity Scale")

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -InfRemoved 161 rows containing missing values (geom_point).Removed 10 rows containing missing values (geom_pointrange).no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -InfRemoved 177 rows containing missing values (geom_point).Removed 10 rows containing missing values (geom_pointrange).Removed 1 rows containing missing values (geom_segment).no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -InfRemoved 519 rows containing missing values (geom_point).Removed 5 rows containing missing values (geom_pointrange).

fig1_row4 <- fig1_fun_v2_row(which_y = "abs_score", ylab = "Absorption")

Removed 15 rows containing missing values (geom_point).Removed 1 rows containing missing values (geom_segment).no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -InfRemoved 269 rows containing missing values (geom_point).Removed 5 rows containing missing values (geom_pointrange).no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -InfRemoved 766 rows containing missing values (geom_point).Removed 5 rows containing missing values (geom_pointrange).

fig1 <- plot_grid(fig1_row1, fig1_row2, fig1_row3, fig1_row4, 
                  ncol = 1)
ggsave("./png_files/fig1.png", plot = fig1, device = "png", width = 12, height = 12 * 0.5)

fig1

Figure 1. Scores on our primary variables of interest—(A) Spiritual Events, (B) Porosity Vignettes, (C) Porosity Scale, and (D) Absorption—for all samples, in all studies. To aid in visual comparison across measures and studies, all scores have been rescaled to range from 0-1. Small points correspond to individual participants, larger points are means, and error bars are ±1 standard deviation; see figure for sample sizes (but note that a few participants in each study were missing data for one or more measures). In Study 1, “faiths of local salience” were as follows—US: Methodism; Ghana: African traditional religion; Thailand: Buddhism; urban China: Buddhism; rural China: spirit mediumship; urban Vanuatu: Presbyterianism; rural Vanuatu: ancestral kastom practices.

Note: the version of this figure included in the main text was edited outside of R: Information about Population (above the plot), Study (below the plot), and sample sizes (in the blank facets) were added manually.

Figure 2 (main text)

fig2_fun <- function(g){
  new_plot <- g +
    geom_point(aes(color = country), alpha = 0.1) +
    geom_smooth(aes(color = country), method = "lm", 
                lty = 2, size = 0.7, alpha = 0, show.legend = F) +
    geom_smooth(method = "lm", color = "black", alpha = 0.7) +
    scale_color_brewer(palette = "Dark2") +
    xlim(0, 1) +
    ylim(0, 1) +
    theme(legend.position = "none") +
    guides(color = guide_legend(override.aes = list(alpha = 1)))
  
  return(new_plot)
}

fig2_study1_pv <- d1_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = pv_score, y = spev_score)) %>%
  fig2_fun() +
  labs(x = "Porosity Vignettes",
       y = "Spiritual Events",
       color = "Country")

fig2_study1_abs <- d1_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = abs_score, y = spev_score)) %>%
  fig2_fun() +
  labs(x = "Absorption",
       y = "Spiritual Events",
       color = "Country")

fig2_study2_por <- d2_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = por_score, y = spev_score)) %>%
  fig2_fun() +
  labs(x = "Porosity Scale",
       y = "Spiritual Events",
       color = "Country")

fig2_study3_abs <- d3_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = abs_score, y = spev_score)) %>%
  fig2_fun() +
  labs(x = "Absorption",
       y = "Spiritual Events",
       color = "Country")

fig2_study4_pv <- d4_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = pv_score, y = spev_score)) %>%
  fig2_fun() +
  labs(x = "Porosity Vignettes",
       y = "Spiritual Events",
       color = "Country")

fig2_study4_por <- d4_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = por_score, y = spev_score)) %>%
  fig2_fun() +
  labs(x = "Porosity Scale",
       y = "Spiritual Events",
       color = "Country")

fig2_study4_abs <- d4_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = abs_score, y = spev_score)) %>%
  fig2_fun() +
  labs(x = "Absorption",
       y = "Spiritual Events",
       color = "Country")

fig2_study1_title <- ggdraw() + 
  draw_label("STUDY 1", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig2_row1 <- plot_grid(
  fig2_study1_title,
  plot_grid(fig2_study1_pv, fig2_study1_abs, ncol = 2, labels = c("A", "B")),
  ncol = 1, rel_heights = c(1, 10))

`geom_smooth()` using formula 'y ~ x'
Removed 20 rows containing non-finite values (stat_smooth).`geom_smooth()` using formula 'y ~ x'
Removed 20 rows containing non-finite values (stat_smooth).Removed 20 rows containing missing values (geom_point).`geom_smooth()` using formula 'y ~ x'
Removed 28 rows containing non-finite values (stat_smooth).`geom_smooth()` using formula 'y ~ x'
Removed 28 rows containing non-finite values (stat_smooth).Removed 28 rows containing missing values (geom_point).

fig2_study2_title <- ggdraw() + 
  draw_label("STUDY 2", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig2_study3_title <- ggdraw() + 
  draw_label("STUDY 3", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig2_row2 <- plot_grid(
  plot_grid(fig2_study2_title, fig2_study3_title),
  plot_grid(fig2_study2_por, fig2_study3_abs, 
            ncol = 2, labels = c("C", "D")),
  ncol = 1, rel_heights = c(1, 10))

`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
Removed 10 rows containing missing values (geom_smooth).

fig2_study4_title <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig2_row3 <- plot_grid(
  fig2_study4_title,
  plot_grid(plot_grid(fig2_study4_pv, fig2_study4_por, 
                      ncol = 1, labels = c("E", "F")), 
            plot_grid(NULL, fig2_study4_abs, NULL, 
                      ncol = 1, rel_heights = c(1, 2, 1), labels = c("", "G", "")), 
            ncol = 2),
  ncol = 1, rel_heights = c(1, 20))

`geom_smooth()` using formula 'y ~ x'
Removed 2 rows containing non-finite values (stat_smooth).`geom_smooth()` using formula 'y ~ x'
Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).`geom_smooth()` using formula 'y ~ x'
Removed 2 rows containing non-finite values (stat_smooth).`geom_smooth()` using formula 'y ~ x'
Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).`geom_smooth()` using formula 'y ~ x'
Removed 2 rows containing non-finite values (stat_smooth).`geom_smooth()` using formula 'y ~ x'
Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 11 rows containing missing values (geom_smooth).

fig_legend <- get_legend(fig2_study1_pv + theme(legend.position = "bottom"))

`geom_smooth()` using formula 'y ~ x'
Removed 20 rows containing non-finite values (stat_smooth).`geom_smooth()` using formula 'y ~ x'
Removed 20 rows containing non-finite values (stat_smooth).Removed 20 rows containing missing values (geom_point).

fig2 <- plot_grid(fig2_row1, fig2_row2, fig2_row3, fig_legend,
                  ncol = 1, rel_heights = c(1, 1, 2, 0.2))

ggsave("./png_files/fig2.png", plot = fig2, device = "png", width = 6, height = 6 * 2.1)

fig2

Figure 2. Relationships between Spiritual Events and measures of porosity (left) and absorption (right), by study and country, rescaled to range from 0-1. Small colored points correspond to individual participants, dotted colored lines correspond to the trend within each country, and solid black lines correspond to the overall trend collapsing across countries. See Figs. S1-S2 for parallel visualizations of other measures of spiritual and secular anomalous events (Studies 3-4).

Figure S1

figs1_study3_abs <- d3_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = abs_score, y = dse_score)) %>%
  fig2_fun() +
  labs(x = "Absorption",
       y = "Daily Spiritual Experience",
       color = "Country")

figs1_study4_pv <- d4_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = pv_score, y = dse_score)) %>%
  fig2_fun() +
  labs(x = "Porosity Vignettes",
       y = "Daily Spiritual Experience",
       color = "Country")

figs1_study4_por <- d4_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = por_score, y = dse_score)) %>%
  fig2_fun() +
  labs(x = "Porosity Scale",
       y = "Daily Spiritual Experience",
       color = "Country")

figs1_study4_abs <- d4_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = abs_score, y = dse_score)) %>%
  fig2_fun() +
  labs(x = "Absorption",
       y = "Daily Spiritual Experience",
       color = "Country")

figs1_study3_title <- ggdraw() + 
  draw_label("STUDY 3", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

figs1_study4_title <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

figs1_row1 <- plot_grid(
  plot_grid(figs1_study4_title, figs1_study3_title),
  plot_grid(figs1_study4_pv, figs1_study3_abs, ncol = 2, labels = c("A", "B")),
  plot_grid(figs1_study4_title, figs1_study4_title),
  plot_grid(figs1_study4_por, figs1_study4_abs, ncol = 2, labels = c("C", "D")),
  ncol = 1, rel_heights = c(1, 10, 1, 10))

`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'

fig_legend <- get_legend(figs1_study4_pv + theme(legend.position = "bottom"))

`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'

figs1 <- plot_grid(figs1_row1, fig_legend,
                   ncol = 1, rel_heights = c(2, 0.2))

ggsave("./png_files/figs1.png", plot = figs1, device = "png", width = 6, height = 6 * 1.2)

figs1

Figure S1. Relationships between Daily Spiritual Experience scores and scores on our measures of porosity (left side) and absorption (right side), by study and country, rescaled to range from 0-1. Small colored points correspond to individual participants (Study 3: N=519; Study 4: N=505), dotted colored lines correspond to the trend within each country, and solid black lines correspond to the overall trend, collapsing across countries.

Figure S2

figs2_study4_pv_hall <- d4_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = pv_score, y = hall_score)) %>%
  fig2_fun() +
  labs(x = "Porosity Vignettes",
       y = "Hallucinations",
       color = "Country")

figs2_study4_por_hall <- d4_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = por_score, y = hall_score)) %>%
  fig2_fun() +
  labs(x = "Porosity Scale",
       y = "Hallucinations",
       color = "Country")

figs2_study4_abs_hall <- d4_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = abs_score, y = hall_score)) %>%
  fig2_fun() +
  labs(x = "Absorption",
       y = "Hallucinations",
       color = "Country")

figs2_study4_pv_para <- d4_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = pv_score, y = para_score)) %>%
  fig2_fun() +
  labs(x = "Porosity Vignettes",
       y = "Paranormal",
       color = "Country")

figs2_study4_por_para <- d4_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = por_score, y = para_score)) %>%
  fig2_fun() +
  labs(x = "Porosity Scale",
       y = "Paranormal",
       color = "Country")

figs2_study4_abs_para <- d4_fig %>%
  # mutate_at(vars(ends_with("_score")), scale) %>%
  ggplot(aes(x = abs_score, y = para_score)) %>%
  fig2_fun() +
  labs(x = "Absorption",
       y = "Paranormal",
       color = "Country")

figs2_study4_title <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

figs2_row1 <- plot_grid(
  plot_grid(figs2_study4_title),
  plot_grid(figs2_study4_pv_hall, figs2_study4_por_hall, figs2_study4_abs_hall,
            ncol = 3, labels = c("A", "B", "C")),
  plot_grid(figs2_study4_pv_para, figs2_study4_por_para, figs2_study4_abs_para,
            ncol = 3, labels = c("D", "E", "F")),
  ncol = 1, rel_heights = c(1, 10, 10))

`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'

fig_legend <- get_legend(figs2_study4_pv_para + theme(legend.position = "bottom"))

`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'

figs2 <- plot_grid(figs2_row1, fig_legend,
                   ncol = 1, rel_heights = c(3, 0.2))

ggsave("./png_files/figs2.png", plot = figs2, device = "png", width = 9, height = 9 * 0.7)

figs2

Figure S2. Relationships between scores on our measures of “secular” anomalous events (A-C: Hallucinations; D-F: Paranormal) and scores on the Porosity Vignettes (left), Porosity Scale (middle), and Absorption (right) scales in Study 4, by country, rescaled to range from 0-1. Small colored points correspond to individual participants (N=505), dotted colored lines correspond to the trend within each country, and solid black lines correspond to the overall trend, collapsing across countries.